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Simulating Heliospheric and Solar Particle Diffusion 
using the Parker Spiral Geometry 

R. C. Tautz, 1 A. Shalchi, 2 and A. Dosch 3 

Abstract. Cosmic Ray transport in curved background magnetic fields is investigated 
using numerical Monte-Carlo simulation techniques. Special emphasis is laid on the So- 
lar system, where the curvature of the magnetic field can be described in terms of the 
Parker spiral. Using such geometries, parallel and perpendicular diffusion coefficients have 
to be re-defined using the arc length of the field lines as the parallel displacement and 
the distance between field lines as the perpendicular displacement. Furthermore, the tur- 
bulent magnetic field is incorporated using a WKB approach for the field strength. Us- 
ing a test-particle simulation, the diffusion coefficients are then calculated by averaging 
over a large number of particles starting at the same radial distance from the Sun and 
over a large number of turbulence realizations, thus enabling one to infer the effects due 
to the curvature of the magnetic fields and associated drift motions. 



1. Introduction 

It is accepted that, in general, the Solar magnetic field can 
be described through an Archimedean spiral as originally 
suggested by Parker [1958]. Today, the Parker model is 
the starting point for countless theoretical — both analytical 
and numerical — investigations. Observations show general 
agreement with the model [Forsyth et al., 1996], especially 
in the ecliptic plane and regarding the most probable field 
direction, although the magnetic field lines are less tightly 
wound as predicted by the model. However, there are overall 
deviations from the Parker structure if the field is measured 
far away from the ecliptic plane, as has been done by the 
Ulysses spacecraft [Smith et al, 2001]. It has been shown 
that such can be explained by the motion executed by the 
foot points of magnetic field lines across coronal loop hole 
boundaries [Schwadron and McComas, 2005], thus giving 
rise to a "sub-Parker" structure. Using numerical simula- 
tions, Riley and Gosling [2007] confirmed the connection of 
a more radial magnetic field and coronal holes. Further- 
more, at high Solar latitudes there are asymmetries in the 
azimuthal field component, and the field lines are generally 
more inclined towards the equator due to the interaction 
with Solar wind plasma [Forsyth et al., 1996]. 

The linear radial decrease of the (azimuthal) field com- 
ponent that was predicted by the Parker model was con- 
firmed from the investigation of Pioneer 10 data [Parker and 
Jokipii, 1976], although some deviations have been known 
depending on the phase of the Solar cycle; for example, 
the radial component is independent on Solar latitude dur- 
ing Solar minima and maxima, and the overall field mag- 
nitude is smaller than expected [Smith, 2004]. Further- 
more, near the Sun the magnetic field has to be described 
through a superposition of dipole, quadrupole, and current 
sheet structures, thus giving rise to extensive modifications 



Centrum fur Physik und Astronomie, Technische 
Universitat Berlin, Hardenbergstrafie 36, D-10623 Berlin, 
Germany. 

2 Department of Physics and Astronomy, University of 
Manitoba, Winnipeg, Manitoba R3T 2N2, Canada. 

3 CSPAR, University of Alabama in Huntsville, Huntsvillc 
AL 35805, USA. 

Copyright 2010 by the American Geophysical Union. 
0148-0227/10/$9.00 



of the Parker structure especially during Solar minima [Ba- 
naszkiewicz et al., 1998]. In defense of Parker, however, 
it should be noted that his model describes the field only 
outside of the zone where field lines execute a rigid-body 
rotation. 

It has been argued that, due to some difficulties with van- 
ishing magnetic field divergence, the standard Parker spiral 
has to be superposed by other field components, preferably 
homogeneous [Goldstein Jr., 1998] with a strength of a frac- 
tion of nT. Such could be interpreted as the local galactic 
magnetic field, for which estimates are available through 
carefully distinguishing the excellently measured Solar field 
and the large-scale field as inferred from pulsar rotation rates 
[Rand and Kulkarm, 1989]. However, it should be noted 
that the form of the magnetic field can be inferred from 
the study of suprathermal electrons (with energies larger 
than 70 eV) , which are aligned with the heliospheric mag- 
netic field (and have therefore been named "strahl"), be- 
cause such a particle distributions broadens with distance 
[Owens et al., 2008]. Independently from other models, such 
leads to a well-formed Parker geometry. 

A modified model [Fisk, 1996] includes the reversion of 
polarity and thus deviates from the Parker model because 
the magnetic field components now depend on the field 
polarity. From the solution of the transport equation, it 
was shown that the energy spectra of protons arriving at 
Earth show significant differences depending on the under- 
lying magnetic field model. Similarly, the transport equation 
has been solved for solar energetic particle events (SEP) us- 
ing both finite difference [Kallenrode, 1993] and stochastic 
methods [Zhang et al. , 2009] . Such particles have sufficiently 
high energy to pose a threat to both astronauts and space 
infrastructure such as satellites. Based on a magnetic field in 
the form of a Parker spiral, the results underline again that 
perpendicular diffusion is important if one wants to under- 
stand the development of the SEP during the time it needs to 
arrive on Earth. Therefore, the knowledge of the magnetic 
field geometry is important if one attempts to understand 
the propagation of particles, which is also important for a 
prediction of "space weather" [Pommois et al., 2002]. 

The smooth, large-scale field described by the Parker — 
and similar, extended — models is, however, only half the 
truth; in reality, the heliospheric magnetic field is highly 
variable on all scales, which can be described in terms 
of waves, fluctuations, or more general turbulence models 
[Smith, 1989]. Thus, the spiral structure becomes evident 
only after (time and space) averaging of the magnetic field, 
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showing that the turbulence strength is of the same order 
of magnitude as the background field. Moreover, as known 
from observations, the small- and middle-scale structures, 
can in turn influence the underlying spiral structure [Dalin 
et al., 2002] by changing the inclination angle of the mag- 
netic field lines. It is generally assumed that the domi- 
nant constituents are magnetohydrodynamic Alfven waves 
[Smith, 1989; Tautz, 2010a], which propagate in opposite 
directions along the background magnetic field lines [Cable 
and Lin, 1998], and convective, two-dimensional structures 
[Schmidt and Marsch, 1995]. However, other waves such 
as magnetosonic waves have also been used, thereby giving 
rise to anisotropies in the power spectra used to characterize 
the turbulence [Chashei, 2000]. Furthermore, magnetosonic 
waves dominate the stochastic acceleration of charged par- 
ticles [Schlickeiser and Miller, 1998; Tautz, 2010a]. 

The motion of charged particles in such systems can 
be described by a diffusion process [see Schlickeiser , 2002; 
Shalchi, 2009, for detailed introductions]. Although the 
problem cannot be considered to be completely solved, sig- 
nificant progress has been made in describing the motion 
of energetic particles such as Cosmic Rays in a turbulent 
medium immersed in a homogeneous background magnetic 
field. For example, simulation results in different turbu- 
lence geometries [Giacalone and Jokipii, 1999; Qin et at, 
2002a, b; Tautz, 2010b] could be reproduced using non- 
linear extensions [e.g., Matthaeus et al., 2003; Shalchi et al., 
2004; Shalchi, 2006; Qin, 2007; Tautz et al, 2008a; Dosch 
et al, 2009; Shalchi, 2010] of the quasi-linear transport the- 
ory [Jokipii, 1966]. 

As soon as it comes to curved background fields, however, 
matters are even less understood. The concept of adiabatic 
focusing has attracted attention recently [e.g., Kunstmann, 
1979; Spangler and Basart, 1981; Bieber et al, 2002; Schlick- 
eiser and Shalchi, 2008], which describes the modification of 
diffusion and transport equations due to a spatial gradient in 
the magnetic field. Other models simply prescribe a curved 
spatial geometry by accepting a background magnetic field 
with non-zero divergence [Minnie et al., 2007]. But a satisfy- 
ing solution has not been achieved yet for real magnetic con- 
figurations that are not uniformly converging or diverging — 
such as the Parker spiral. In principle, the real field lines, 
which are distorted by the presence of turbulence, can be 
inferred using a field line wandering (or field line random 
walk) approach [Jokipii and Parker, 1969; Matthaeus et al., 
1995; Shalchi and Kourakis, 2007]. By specifying the turbu- 
lence power spectrum, a direct calculation of the stochastic 
field lines is then possible, showing excellent agreement with 
Helios 2 data. Furthermore, Kobylmski [2001] showed that 
the diffusion tensor can be calculated analytically by apply- 
ing two rotations to the parallel and perpendicular diffusion 
coefficients, although, in general, such a transformation will 
always be ambiguous. 

Furthermore, the question of "diffusivity" , i.e., whether 
the diffusion coefficients attain finite values, is still not 
solved satisfactorily. Whereas, in perpendicular scatter- 
ing in magnetostatic slab turbulence (in Alfvenic slab tur- 
bulence, diffusion is recovered Shalchi et al. 2007; Tautz 
2010a), sub-diffusion is clearly established, it is generally 
assumed that, in isotropic turbulence, the particle motion 
is diffusive — although there are contradicting analytical re- 
sults [Tautz et al, 2008b]. However, for anisotropic non- 
slab turbulence, such is not clear [Zimbardo et al., 2006]. 
There are several cases where sub- or even super-diffusion 
has been found in numerical simulations [Tautz and Shalchi, 
2010] and in observations of electrons accelerated at inter- 
planetary shocks [Perri and Zimbardo, 2007]. It has even 
been thought that the underlying process might not be a 
classic diffusion process but a Levy random walk [Zimbardo, 
2005]. However, such analyses would first require a stable 
description of the turbulence, because the results depend 



sensitively on the details of the turbulence model and its 
parameters. 

In this article, the transport of test-particles is investi- 
gated using a numerical Monte Carlo technique [Giacalone 
and Jokipii, 1999; Michalek, 2001; Tautz, 2010b]. The mag- 
netic field is composed of a large-scale Parker-type field and 
a turbulent component, which is assumed to be isotropic. It 
is known from observations [Matthaeus et al, 1990; Bieber 
et al., 1994; Chashei, 2000] that the turbulent magnetic field 
component is most likely not isotropic. There are even stud- 
ies that incorporate such effects in sophisticated turbulence 
models [Sridhar and Goldreich, 1994; Goldreich and Srid- 
har, 1995]. The problem, however, is that the basic trans- 
port theories are not sufficiently understood yet to allow for 
the isolation of small effects. Here, therefore, an isotropic 
turbulence spectrum will be used [Tautz et al., 2006, 2008a] 
that uses a constant energy range. Numerically, the well- 
known method will be applied of superposing a homogeneous 
(in our case: the Parker spiral) magnetic field with a tur- 
bulent component calculated by the summation over plane 
waves with random directions of propagation and random 
phase angles [Tautz, 2010b]. Thereby, the diffusion coef- 
ficients are calculated from the mean square deviation. It 
will be shown how "parallel" and "perpendicular" particle 
displacements must be transformed using concepts from dif- 
ferential geometry. The numerical ansatz will be relativisti- 
cally correct, although such might not be necessary for the 
energies considered. In a second paper [Tautz and Vocks, 
2010], the effects of Whistler wave turbulence will be stud- 
ied, which is important for the formation of the "strahl" elec- 
trons that have been mentioned above. Here, in contrast, 
magnetostatic turbulence will be employed. Although, as 
has already been mentioned, plasma waves may have a se- 
vere impact on the diffusion coefficients especially if their 
electric fields are of significant magnitude, such a configura- 
tion represents the most basic case that has to be understood 
first before implementing more advanced features. The ap- 
proach is, to some extent, comparable to a method used by 
Pei et al. [2006], although in their work focus was laid on 
the calculation of SEP onset times. It was found that, in 
the presence of turbulent field component, the onset times 
were generally reduced in comparison to a smooth Parker 
magnetic field. 

The present paper is organized as follows: In Sec. 2, 
the geometry of the Parker spiral is introduced and it is 
shown how appropriate modifications can be made to the 
(numerical) implementation of the simulation code that cal- 
culates the diffusion coefficients. In Sec. 3, the setup of 
the test-particle simulations and the turbulence generation 
are briefly explained. In Sec. 4, the simulation results are 
presented and compared both to previous calculations and 
simulations in homogeneous background magnetic fields. Fi- 
nally, Sec. 5 provides a summary and conclusions regarding 
future work. 

2. Parker Spiral Geometry 

According to Parker [1958] , the large-scale spiral pattern 
of the Solar magnetic field in interplanetary space (without 
taking into account the turbulent component) is determined 
through an equation relating the radial distance from the 
Sun, r, and the azimuth angle in the ecliptic plane, <f>, as 

where the parameter ip is the azimuth angle of the field line 
at the co-rotation radius. Furthermore, vsw ~ 400km/s is 
the outward velocity of the solar wind, ujq is the angular 



TAUTZ ET AL.: PARKER SPIRAL GEOMETRY 



X- 3 



velocity of the sun, and b = 46 x 10 AU is a distance be- 
yond which any direct influence of the sun may be neglected 
[cf. Parker, 1958]. Physically, 6 marks the end of the region, 
where the field lines execute a "rigid rotation" and, there- 
fore, point radially outwards; i. .e., for r < b the field lines 
point radially outward while they corotate with the Sun's 
surface and, thus, Eq. (1) is not valid. Here, however, only 
the case r > b will be considered. 

In what follows, quantities and relations will be referred 
to Figs. 1 and 2, which illustrate the explanations given 
here. Solving Eq. (1) for the azimuth angle as a function 
of the radius, r, yields 



4>{r) 



b rr 
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with <" = vsw/uq ~ 1 AU [Zank et al, 2004]. 

Let us now consider an arbitrary particle in the ecliptic 
plane (black dot) at the cartesian coordinates (a; s , y s ) denot- 
ing the starting point on a magnetic field line (black dash 
dotted line). The starting point may be expressed through 
the new coordinates (r s ,<fi s ) using the transformations 



r s = V x 1 + y'i 



= arctan 



(3a) 
(3b) 



where, on evaluating the arctan function, attention has to 
be paid to choosing the right quadrant. The inner field line 
azimuth angle, ip, can be determined by back-tracing the 
field line to the radius r = b. Since the starting point lies 
on that field line, ip can simply be calculated by 



ip = <p s 



In 



(4) 



which is visualized in Fig. 1. 

Usually the end point of the particle trajectory, (x e ,y e ), 
will be located on a different (but nearby) field line (see 
brown dot on the brown dashed line in Fig. 2). Particu- 
larly, the particle has drifted from one field line to another 
due to diffusion processes. According to the definition of the 
parallel and perpendicular (running) diffusion coefficients, 
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one has to calculate the parallel and perpendicular mean 
square displacements, respectively, which are denoted by 
{A 2 }. The perpendicular displacement is associated with the 
shortest distance — perpendicular to the field lines — between 
the ending point and the initial field line, leading to the tar- 
get point (red dashed line and red dot). In a uniform mag- 
netic field such is simply a straight line perpendicular to 
the background magnetic field lines. Likewise, the parallel 
displacement is the shortest distance between the starting 
point and the target point; in a uniform magnetic field, it 
would be a straight line, too. But how does that have to be 
translated for a curved background magnetic field such as 
that described by the Parker spiral? 

The natural answer is motivated by differential geometry. 
The perpendicular displacement is characterized by a curve 
perpendicular to all field lines, which runs through the end 
point and the target point. As a first approach it is as- 
sumed that such a curve can be approximated by a straight 
line as in Euclidean geometry (red dashed line), based on the 
fact that the turbulence in interplanetary space is moderate 
[Bieber et al, 1994] and, therefore, perpendicular diffusion 
is weak, too. Furthermore, it was shown by Shalchi and 
Dosch [2009] that, even in strong isotropic turbulence, one 
generally has k± <C k\\ . 

This also motivates the assumption that the particle tra- 
jectory ends on nearby field lines and states that the parallel 



displacement is given through the arc length of the initial 
(black dash dotted line) field line to a (red dot) target point, 
where the shortest distance to the (brown dot) particle end 
point is perpendicular on all field lines. Hence, instead of a 
straight line the "shortest distance" is now a geodesic fol- 
lowing the geometry of the Parker spiral. As an approxi- 
mation, however — owed to the fact that the turbulence in 
interplanetary space is weak and that, therefore, perpen- 
dicular diffusion is weak, too — the calculation of the target 
point proceeds as follows: One simply calculates the point 
from which the distance — this time taken as a straight line 
as in Euclidean geometry — to the particle trajectory's end 
point is minimal (red dashed line). Then, the perpendicular 
displacement, denoted by As±, is simply given through the 
distance between the target point, r t , and the end point, r e , 
as 

As\= min I [r t cos(0(r t )) - x e ] 2 + [r t sin(0(n)) - y e ] 2 \, 

(6) 

which function is visualized in Fig. 3. 

The only intricacy is to find the correct minimum, since 
the function in Eq. (6) has multiple local minima for <f>+2irk 
with k approximately an integer number. Hence, a method 
in three steps will be used. First step: to identify a minimum 
of any function /, three points a < b < c are needed with 
f(b) < f(a) and f(b) < /(c) are necessary — which is called 
"bracketing" of the minimum [Press et ai, 2007, p. 490]. 
Assuming p t ~ <p c and assuming that the particle will not 
complete one or more full orbits around the Sun, one can 
therefore search for a minimum starting from r B . By itera- 
tively solving for r{<p e ~ <pt) starting with ro = r B one has, 
therefore, 

n+i=fe[l + ln(^)] +C(0c-V), (7) 

which will be stopped after a few steps. For very small dis- 
placements (r c ,(f> c ) ~ (r s ,<fis), i.e., after short times, how- 
ever, the iteration can yield wrong results. Second step: In 
that fail-safe method will determine an approximate 

bracketing condition from calculating r(<p B — 2ir) < r B < 
r((f> B + 2tt), where again use is made of the assumption that 
a given particle will not complete a full orbit around the 
Sun. Third step: Based on the bracketing of the minimum, 
it is then a simple task to calculate the exact minimum us- 
ing, e.g., Brent's method [Press et at, 2007, pp. 496-502] 
without derivative information [see also Forsythe et al. 1977, 
§8.2, Brent 2002, Chpt. 5]. 

Once the target point has been determined, the parallel 
displacement, As\\, is calculated through the arc length of 
the field line between the target point and the starting point, 
yielding 



As H =C( 



C(r B 



(8) 



where the arc length function, C, can be calculated analyt- 
ically and is given through 



C { r) = [ d sJl + s* 



r - b y/(r - b) 2 + C 2 C A . ,/r-6 

~y~ c ~ 2 Arsmh [~r 



(9a) 
(9b) 



For instance, the trajectory parameters from Fig. 2 (with 
( = l)have been used to obtain As\\/b = 20.2 [by using 
Eq. (8)] and As ± /b = 2.2. 

Until here, only the ecliptic plane — which, in the simu- 
lation box, corresponds to the x-y plane — has been consid- 
ered. To incorporate drift motion, the vertical displacement, 
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Az = z c — z s , is calculated so that, from the comparison of 
As± and A z , the drift effect can be inferred. 

Furthermore, the (background) magnetic field compo- 
nents, which are designed to cover mainly the ecliptic 
plane, are described using spherical coordinates through [see 
Parker, 1958; Burger and Hitge, 2004] 



B [ b - 

r 



B r (r,e,4 
B e {r,e,cf>) = 



B.Xr.O.o) = a, ( ^ r —^- sinf 



(10a) 
(10b) 
(10c) 



with Bo = 1830 nT [Zank et al, 2004]. Note that the 
strength of the magnetic field is maximal in the ecliptic 
plane, and that no care has been taken of the current sheet 
and the polarity reversal [cf. Fisk, 1996; Burger and Hitge, 
2004]. For small r, the field is purely radial and scales with 
r~ , whereas, according to Eqs. (10), for large r the field 
becomes more and more azimuthal and scales with in 
agreement with Pioneer 10 data [Parker and Jokipii, 1976]; 
see also Fig. 2 and Fig. 6 of Parker [1958]. 

When the diffusion coefficient is calculated in the way 
described here, the curved space denned by the large-scale 
magnetic field is mapped to a Euclidean space. The diffu- 
sion coefficient can include the additional terms coming out 
of the divergence in the curved space. Thus, the diffusion 
coefficient carries a meaning different from local diffusion 
coefficient defined in the curved space. Note also that, until 
here, only the smooth spiral field lines have been considered. 
Both in the real world and in the simulation code, the back- 
ground field will be superposed by a turbulent component. 
However, diffusion coefficients will be calculated in reference 
to the mean field component. 



3. Monte-Carlo Simulations 

For the simulation of cosmic ray scattering processes, 
a modified version of the recently developed Padian code 
[Tautz, 2010b] has been used, which traces the trajectories 
of a large number of test particles (typically 10 4 ) for a suf- 
ficiently long time (typically 10 4 Larmor orbits). In doing 
so, it can be decided whether the transport is diffusive, un- 
perturbed, or subdiffusive, and the diffusion coefficients can 
be determined. For the integration of the equation of mo- 
tion, an integration algorithm with adaptive step sizes such 
as the Bulirsch-Stoer method [see Stoer and Bulirsch, 2002; 
Press et al., 2007, pp. 921-928] is used, which limits the rel- 
ative deviation in the particle rigidity, AR/R, to be smaller 
than 10 _3 %. Here, the particle is expressed as momentum 
per unit charge times speed of light, i.e., R = pc/e. Note 
the difference to several previous articles [e.g., Tautz et al., 
2008a; Tautz, 2010b], where a dimensionless "rigidity" -like 
quantity was defined through R = 'ymv/((2£o), where fi 
denotes the gyro-frequency. 

In comparison to the original formulation [Tautz, 
2010b], the equation of motion — i.e., the Newton-Lorentz 
equation — has to be slightly modified due to the non- 
constant background magnetic field. Furthermore, SI units 
are used so that the rigidity and the magnetic field strengths 
can be given in mega Volt (MV) and nanoTesla (nT), respec- 
tively. The trajectory x(t) (normalized to a characteristic 
length scale £o; see below) is therefore calculated as a func- 
tion of the dimensionless time r = vt/£o using 



d x 

d7 T 



R 



R = |; R x [B (r) e Bo + SB(r) e 6B ] , 



(11a) 
(lib) 



where Eq. (lib) corresponds to p — qv x B (SI units). The 
scaling factor is given by a = £oc/10 15 ~ 1345.45 m /s for 
typical parameters (see below). Note that Eqs. (11) are 
valid for all particle species and that the rest frame of the 
Sun has been used. 

The unit vectors es and ess point in the directions of 
the Parker field and the turbulent magnetic field, respec- 
tively. In the ecliptic plane, the total background mag- 
netic field strength is, approximately, determined through 
Eq. (10) as 



B (r) = 1830 nT 



(12) 



Following Zank et al. [1996, Sec. 3.1], a WKB approxima- 
tion [see also Matthaeus et al, 1994] leads to an equation 
for the turbulent magnetic field energy per mass, Esb = 
8B 2 /(8irp), as 

^ + ^ = 0, (13) 
dr r 

yielding Esb oc r _1 . Equation (13) was motivated by the 
Walen relation [Matsuoka et al, 2002], i.e., AB oc pAv, 
where AB and Av are changes in the magnetic field and 
the plasma velocity, respectively. Based on the Walen re- 
lation, the ratio of kinetic and magnetic energy (both per 
mass) are constant. Therefore, with the additional assump- 
tion that, in the ecliptic plane, the mass density behaves as 
par" 



SB 

8B rc 



(14) 



where 5B TC f — 4nT is the turbulent field strength at r re f = 
1 AU. Accordingly, the turbulent magnetic field strength is 
given through [Zank et al, 2004; Shalchi et al., 2010, Ap- 
pendix B] 

8B(r) = 4nT(r[AU]r 3/2 , (15) 

where r is the radial distance from the Sun in astronomi- 
cal units (AU). Based on Eq. (12), the relative strength of 
turbulent and background magnetic field depends on r as 



SB_ 

Bo 



'[AU] 



1 + (r [AU])" 



(16) 



The turbulence model, which, to some extent, is based 
on a model used by Giacalone and Jokipii [1999], gener- 
ates random magnetic fluctuations by the superposition of a 
large number N of plane waves (typically 512) with random 
directions of propagation, random phases, and amplitudes 
as prescribed by the turbulence spectrum of Tautz et al. 
[2006, 2008a] that has the form 



G(k) = [1 + (lokf 



iok) 



-5/3 



(17) 



where the parameter v = 5/6 denotes the Kolmogorov spec- 
tral index for a power-law k~ 5 ' 3 [cf. Podesta et al, 2007] 
and where to is the isotropic turbulence bend-over scale. 
Note that, in the Solar system, the bend-over scale, to, may 
depend on the position [Bruno and Carbone, 2005]. The 
spectrum from Eq. (17) has a constant energy range (i.e., 
the range k < tg 1 ). For a more general spectrum, a formula 
similar to that of Eq. (17) can be found, e.g., in Shalchi 
and Weinhorst [2009]. 

The spatially fluctuating but time-independent (i.e., 
magnetostatic) magnetic field is then calculated as 



5B(x, y, z) = ReJ2 e' ± ^(fen)e i ( fe " 



:'+Pn) 



(18) 
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where the sum extends over N logarithmically spaced 
wavenumber values k n . The amplitude function A(k n ) oc 
sj G(k n ) is related to the turbulence spectrum and e' ± is 
a unit vector in the direction perpendicular to z . The 
primed coordinates are obtained from a rotation matrix, 
whose angles are randomly generated for each summand n. 
The parameter /3 is the random phase for each wave mode. 
Thus, Eq. (18) generates an isotropic, time-independent 
(i.e., magnetostatic) turbulent magnetic field. Accordingly, 
the orientation of (x,y,z) is irrelevant; however, for conve- 
nience the z direction is kept perpendicular to the ecliptic 
plane. The divergence of SB is kept zero due to the fact that 
e'j_ _L fe n for all summands, which corresponds to V -SB = 
in Fourier space. 

Care must be taken about the minimum and maximum 
wave numbers. Whereas the analytical form of the spec- 
trum from Eq. (17) extends over all wavenumbers, such is 
not possible in computer simulations. There are two major 
conditions that have to be fulfilled, which are: (i) the reso- 
nance condition stating that there must exist a wavenumber 
so that kjj,Rh = 1 is fulfilled; and (ii) the time scale condi- 
tion stating that J?t m axfcmm-RL < 1- The second condition 
requires a maximum turbulence scale (defined through 
to be larger than the distance traveled by the particle to en- 
sure the particle cannot move out of the system. 

From simulations in a homogeneous background magnetic 
field [e. g., Tautz, 2009], one knows that both conditions (i) 
and (ii) do not depend on the turbulence bend-over scale 
£o; instead, they depend on the maximum and minimum 
scale of the system (given by the minimum and maximum 
wavenumber, respectively). Condition (ii) therefore states 
that the particles must not travel farther than the system 
scale L max = fc~f n . Although, in the simulation code, the 
turbulence is generated wherever the particle position is, 
particles start to free-stream once the condition is violated. 
In that case, one finds k± oc t -1 , which is equivalent to a 
constant perpendicular mean square displacement and indi- 
cates a purely parallel motion. 

In Fig. 4, the magnetic field strength is shown as a func- 
tion of the distance from the Sun. In the Padian simula- 
tion code, the field strength is usually normalized to unity; 
here, however, the field strength decreases with increasing 
distance. Furthermore, the (average) magnitude of the tur- 
bulent field is shown in comparison to the background field 
strength, indicating a maximum at 1 AU. 

By integrating the particle trajectory as given by 
Eq. (11), the mean square displacement is then determined 
as described in Sec. 2. Finally, the diffusion coefficients 
and the mean free paths are calculated by averaging over all 
particles and all different turbulence realizations. 

4. Simulation Results and Comparison 

For the simulation runs, the particle energies were cho- 
sen in the range 10 1 MV < R < 10 4 MV. All runs were 
carried out using a total number of 1000 particles in 25 
different turbulence realizations. It is interesting that, in 
contrast to simulations involving homogeneous background 
magnetic fields, here small particle energies require consid- 
erably longer computation times than high energies. 

Furthermore, all particles are assumed to start in the 
ecliptic plane with an equal initial radial distance to the 
sun of r s = 10£o- Since the bend-over scale is usually taken 
to be £o — 0.03 AU, a heliocentric distance of 1AU corre- 
sponds to 33£o. While, in the absence of a turbulent mag- 
netic field component, the particles follow the magnetic field 
lines, there is extensive scattering if the magnetic field is tur- 
bulent (see Fig. 5). Especially at low particle energies, there 
even seems to be confinement-like mechanisms (not shown in 



the figure), where particles repeatedly move back and forth 
along the same magnetic field line without covering a large 
net distance. 

In Fig. 6, the resulting mean free paths in the directions 
parallel, perpendicular, and vertical to the Sun's magnetic 
field are shown as a function of the normalized dimension- 
less time, vt/£o- (Note that, here, the term "perpendicular" 
is used to describe a direction in the ecliptic plane, whereas 
"vertical" is oriented normal to the ecliptic plane.) 

The comparison of the perpendicular and vertical mean 
free paths shows no qualitative difference, which leads to the 
conclusion that a drift effect can be neglected. The result 
that turbulent magnetic fields suppress particle drift mo- 
tions is supported by the work of Minnie et al. [2007], who 
found that, if scattering processes are important, all drift 
motions are superseded by diffusion. 

Furthermore, note that all three mean free paths in Fig. 6 
seem to approach a constant, thus indicating a diffusive 
behavior. In contrast to a configuration with a homoge- 
neous background magnetic field as investigated in Tautz 
and Shalchi [2010], there is no indication here of a (even 
slightly) subdiffusive behavior — at least not for sufficiently 
high rigidities R > 10 MV. 

In Fig. 7, the parallel mean free path is shown as a 
function of the particle rigidity. Furthermore, a compari- 
son is shown to the Palmer consensus range (shaded box) 
consisting of mean free path values from various observa- 
tional studies [see Palmer, 1982; Bieber et al, 1994, and 
references therein]. 

Note that all simulation results were obtained for an ini- 
tial radial distance r B = 10£o- However, the difference to 
additional simulation runs with r s — 30£o turned out to be 
negligible. Therefore, the radial dependence of the magnetic 
field strength seems to have only a small effect on the mean 
free paths — in contrast to the turbulent magnetic field. Ac- 
cording to quasi- linear theory [e.g., Jokipii, 1966], the tur- 
bulence strength relates to the parallel mean free path as 
Bo/\8B\ oc v/Ay", which behavior has also been observed 
in simulations [Tautz, 2010a] for intermediate turbulence 
strengths; for strong turbulence, such may be different [cf. 
Shalchi et al., 2009]. 

In Fig. 8, the ratio of the perpendicular/vertical and par- 
allel mean free path is shown as a function of the particle 
rigidity. In agreement with previous work [e.g., Giacalone 
and Jokipii, 1999; Shalchi et al., 2004; Zank et al., 2004; 
Tautz et al., 2006], Ax/Ay tends to be slightly reduced as the 
particle rigidity increases. The overall values are somewhat 
higher than normally observed, because, using the estimate 
of Bieber et al. [2004], i.e., Ax/Ay « (SB/B ) 2 /6, one has, at 
a radial distance of 10 and 20 AU yields values for A/x/Ay of 
approximately 0.05 and 0.07, respectively, whereas the sim- 
ulation results are all well above Ax/Ay > 0.1 (see Fig. 8). 

Finally, the ratio of vertical and perpendicular mean free 
paths is shown in Fig. 9, illustrating that, although there is 
are noticeable uncertainties, the ratio is in agreement with 
unity. Therefore, the drift effects can be considered as neg- 
ligible as already pointed out by Minnie et al. [2007]. 

5. Summary and Conclusion 

In this article, particle scattering and diffusion param- 
eters have been calculated for the Solar system. The in- 
vestigation was based on the Parker spiral model for the 
background magnetic field together with a WKB model to 
account for the turbulent magnetic field component, which 
varies with the radial distance. For the first time, the Parker 
geometry of the Sun's magnetic field has been implemented 
in a test-particle simulation code. 

The prediction of negligible drift motions could be con- 
firmed (Fig. 9). It is obvious from Fig. 7, the results for the 
parallel mean free path shown here are too small, except for 
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high rigidities R > 10 4 MV. There are three possible expla- 
nations: 

First, in the present version an isotropic turbulence model 
has been used. Usually, a slab/2D model consisting of two 
turbulence components with wavevectors parallel and per- 
pendicular to the background magnetic field are used for 
applications in the Solar system [Bieber et al, 1994, 1996], 
which increases the parallel mean free path [e.g., Tautz 
and Shalchi, 2010]. On the other hand, it has been shown 
by Weinhorst and Shalchi [2010] that a more realistic tur- 
bulence model, which overcomes the singularities inherent 
in the slab/2D formulation, can reproduce data calculated 
from spacecraft measurements better. However, the turbu- 
lence model that was designed in the style of the so-called 
Maltese cross [Matthaeus et al, 1990] has not yet been im- 
plemented into simulation codes. 

Second, the WKB model that has been used for the 
strength of the turbulent magnetic field could overestimate 
the true turbulence level. To allow for better adjustment of 
the turbulence parameters, more measurements of the tur- 
bulence level at different radial positions would be necessary. 

Third, for other choices of the correlation scales, espe- 
cially for different scales in the directions parallel and per- 
pendicular to the background magnetic field, the resulting 
transport parameters would be different [see also Pet et al., 
2010]. Especially, a different choice of parallel and perpen- 
dicular turbulence bend-over scales could modify the result- 
ing mean free paths [Matthaeus et al., 2003; Shalchi, 2009, 
Sec. 5.4.4]. 

Future work should, therefore, try to implement a Mal- 
tese cross based turbulence model to account for the various 
anisotropies between the parallel and the perpendicular tur- 
bulent magnetic field components [e.g., Narita and Glass- 
meter, 2010; Chen et at, 2010]. The effects due to the Solar 
wind would also be worth studying. Furthermore, a reli- 
able description is necessary of the turbulent magnetic field 
strength on the radial distance from the Sun. 
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Figure 1. Illustration of the angles used to describe the 
Parker spiral. The parameter tp corresponds to the an- 
gle at which the large-scale field line (solid line) starts, 
whereas the variable (f>(r) determines the course of the 
field line in dependence on the radius, r. 
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Figure 2. (Color online) The meaning of "parallel" 
and "perpendicular" in Parker spiral geometry. Start- 
ing (black dot) on the black dash dotted field line, the 
particle trajectory (blue solid line) ends (brown dot) on 
the brown dashed field line. As an approximation, the 
perpendicular distance between the two field lines (red 
short-dashed line) is taken to be a straight line, where 
the target point (red dot) has to be calculated through 
minimization of the perpendicular distance. 
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Figure 3. The function from Eq. (6), the minimum 
of which (black dot) has to be found in order to calcu- 
late the radius of the target point, r t , and the associated 
perpendicular and parallel displacements, As± and Aa\\, 
respectively. The correct minimum is that next to the 
starting point, for which r E /b — 4.47. 
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Figure 4. The decrease of the background magnetic 
field strength with radial distance from the Sun as pre- 
scribed by the Parker spiral through Eqs. (10). At 
r — r s = 33^o, one has Bo = 1.68 nT. The inset shows 
the relative strength of the turbulence in terms of the 
background magnetic field. 
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Figure 5. Particle trajectories without (upper panel, 
dashed lines) and with (lower panel, solid lines) magnetic 
turbulence. As a guide to the eye, the magnetic field lines 
according to the Parker spiral (with, in the lower panel, 
superposed by the turbulent component) are shown as 
dotted lines. Furthermore, the dash dotted line roughly 
marks the Earth orbit around the Sun. 
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Figure 7. The mean free path as a function of the par- 
ticle rigidity. The error bars indicate the simulation re- 
sults, and the shaded area illustrates the Palmer [1982, 
see also Bieber et al. 1994] consensus range. 
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Figure 8. The ratio of the perpendicular and the parallel 
mean free path (upper panel) and the ratio of the vertical 
and the parallel mean free path (lower panel). The error 
bars are calculated from the individual estimated mean 
errors of the individual mean free path components. 
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Figure 6. The mean free paths normalized to the astro- 
nomical unit, k/(v£o), as a function of the dimensionless 
simulation time, vt/lo. The solid line shows the paral- 
lel mean free path (i.e., the coefficient of the diffusion 
along the curved magnetic field lines of the Parker spi- 
ral), whereas the dashed line show the mean free path 
perpendicular but in the ecliptic plane. The dotted 
line shows the vertical mean free path normal to the 
ecliptic plane, which is defined in the classic way, i. e., 
\ z = 3{(Az) 2 )/(2vt). 
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Figure 9. The ratio of the perpendicular and the parallel 
mean free path (upper panel) and the ratio of the vertical 
and the parallel mean free path (lower panel). The error 
bars are calculated from the individual estimated mean 
errors of the individual mean free path components. 



